Making maps with R

Instituto de Medicina Molecular & Immune-Stats

Description

Libraries

library(here)
library(data.table)
library(magrittr)
library(ggplot2)
library(scales)
library(colorspace)
library(janitor)
theme_set(theme_bw())

1 Making maps with R

• link

Data from the North Carolina data

library(sf)
nameshp <- system.file("shape/nc.shp", package = "sf")
d <- st_read(nameshp, quiet = TRUE)
# transform
d <- st_as_sf(d)

# infant deaths 1974
d$vble <- d$SID74
# infant deaths 1979
d$vble2 <- d$SID79
head(d)

## Simple feature collection with 6 features and 16 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: -81.74107 ymin: 36.07282 xmax: -75.77316 ymax: 36.58965
## Geodetic CRS:  NAD27
##    AREA PERIMETER CNTY_ CNTY_ID        NAME  FIPS FIPSNO CRESS_ID BIR74 SID74
## 1 0.114     1.442  1825    1825        Ashe 37009  37009        5  1091     1
## 2 0.061     1.231  1827    1827   Alleghany 37005  37005        3   487     0
## 3 0.143     1.630  1828    1828       Surry 37171  37171       86  3188     5
## 4 0.070     2.968  1831    1831   Currituck 37053  37053       27   508     1
## 5 0.153     2.206  1832    1832 Northampton 37131  37131       66  1421     9
## 6 0.097     1.670  1833    1833    Hertford 37091  37091       46  1452     7
##   NWBIR74 BIR79 SID79 NWBIR79                       geometry vble vble2
## 1      10  1364     0      19 MULTIPOLYGON (((-81.47276 3...    1     0
## 2      10   542     3      12 MULTIPOLYGON (((-81.23989 3...    0     3
## 3     208  3616     6     260 MULTIPOLYGON (((-80.45634 3...    5     6
## 4     123   830     2     145 MULTIPOLYGON (((-76.00897 3...    1     2
## 5    1066  1606     3    1197 MULTIPOLYGON (((-77.21767 3...    9     3
## 6     954  1838     5    1237 MULTIPOLYGON (((-76.74506 3...    7     5

2 ggplot2

library(ggplot2)
library(viridis)
# ggplot(d) +
#   geom_sf(aes(fill = vble)) +
#   scale_fill_viridis() + 
#   theme_bw()
ggplot(d) +
  geom_sf(aes(fill = vble)) +
  scale_fill_continuous_sequential() + 
  theme_bw() +
  theme(
    legend.position = "top"
  )

2.1 plotly

library(plotly)
g <- ggplot(d) +
  geom_sf(aes(fill = vble), color = "gray5", size = 0.1) +
  scale_fill_continuous_sequential() + 
  theme_bw()
ggplotly(g)

3 leaflet

The sf object that we pass to leaflet() needs to have a geographic coordinate reference system (CRS) (EPSG code 4326) indicating latitude and longitude. Here, we use the st_transform() function of sf to transform the data d which has CRS given by EPSG code 4267 to CRS with EPSG code 4326.

st_crs(d)

## Coordinate Reference System:
##   User input: NAD27 
##   wkt:
## GEOGCRS["NAD27",
##     DATUM["North American Datum 1927",
##         ELLIPSOID["Clarke 1866",6378206.4,294.978698213898,
##             LENGTHUNIT["metre",1]]],
##     PRIMEM["Greenwich",0,
##         ANGLEUNIT["degree",0.0174532925199433]],
##     CS[ellipsoidal,2],
##         AXIS["latitude",north,
##             ORDER[1],
##             ANGLEUNIT["degree",0.0174532925199433]],
##         AXIS["longitude",east,
##             ORDER[2],
##             ANGLEUNIT["degree",0.0174532925199433]],
##     ID["EPSG",4267]]

library(leaflet)

pal <- colorNumeric(palette = "YlOrRd", domain = d$vble)

l <- leaflet(d) %>% 
  addTiles() %>%
  addPolygons(
    color = "white",
    fillColor = ~ pal(vble),
    fillOpacity = 0.8
  ) %>%
  addLegend(pal = pal,
            values = ~ vble,
            opacity = 0.8)
l

l %>% addMiniMap()

3.1 Save interactive figure as png

# Saves map.html

library(htmlwidgets)
saveWidget(widget = l, file = here::here("results/figures", paste(Sys.Date(), "map.html", sep = "_")))

# Takes a screenshot of the map.html created and saves it as map.png
library(webshot)
# webshot::install_phantomjs()
webshot(url = here::here("results/figures", paste(Sys.Date(), "map.html", sep = "_")), file = here::here("results/figures", paste(Sys.Date(), "map.png", sep = "_")))

4 mapview

library(mapview)
mapview(d, zcol = "vble")

library(RColorBrewer)
pal <- colorRampPalette(brewer.pal(9, "YlOrRd"))
# version 2
mapview(d,
        zcol = "vble",
        map.types = "CartoDB.DarkMatter",
        col.regions = pal)

# version 3
map1 <- mapview(d, zcol = "vble")
leaflet::addMiniMap(map1@map)

5 Side by side plots with mapview

library(leaflet.extras2) # similar to patchwork
m1 <- mapview(d, zcol = "vble")
m2 <- mapview(d, zcol = "vble2")
m1 | m2

6 Synchronized maps with leafsync

m1 <- mapview(d, zcol = "vble")
m2 <- mapview(d, zcol = "vble2")
m <- leafsync::sync(m1, m2)
m

htmltools::save_html(html = m, 
                     file = here::here("results/figures", paste(Sys.Date(), "map-side-by-side.html", sep = "_")))

7 tmap

library(tmap)

# define the mode
tmap_mode("plot") # static
# tmap_mode("view") # interactive

tm_shape(d) + 
  tm_polygons("vble") #+

  # tm_shape(st_centroid(d)) + 
  # tm_dots("vble")

8 Mobility flows with flowmapblue

• To map mobility data
• data from cellphone companies

# devtools::install_github("FlowmapBlue/flowmapblue.R")
library(flowmapblue)

# data frame 1
locations <- data.frame(
  id = c(1, 2, 3),
  name = c("New York", "London", "Rio de Janeiro"),
  lat = c(40.713543, 51.507425, -22.906241),
  lon = c(-74.011219, -0.127738, -43.180244)
)
locations

##   id           name       lat        lon
## 1  1       New York  40.71354 -74.011219
## 2  2         London  51.50742  -0.127738
## 3  3 Rio de Janeiro -22.90624 -43.180244

# atadate frame 2
flows <- data.frame(
  origin = c(1, 2, 3, 2, 1, 3),
  dest = c(2, 1, 1, 3, 3 , 2),
  count = c(42, 51, 50, 40, 22, 42)
)
flows

##   origin dest count
## 1      1    2    42
## 2      2    1    51
## 3      3    1    50
## 4      2    3    40
## 5      1    3    22
## 6      3    2    42

# final plot of transitions fromand to the 3 destinations
# flowmapblue(locations = locations, flows = flows, mapboxAccessToken = NULL, clustering = TRUE, darkMode = TRUE, animation = TRUE)